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ABSTRACT 

We investigate both analytically and numerically the motion of massless particles orbiting primary star in a close circular binary 
system with particular focus on the gas drag effects. These are the first calculations with particles ranging in size from 1 m to 10km, 
which account for the presence of a tidally perturbed gaseous disk. We have found numerically that the radial mass transport by the 
tidal waves plays a crucial role in the orbital evolution of particles. In the outer region of the gaseous disk, where its perturbation is 
strongest, the migration rate of a particle for all considered sizes is enhanced by a factor of 3 with respect to the axisymmetric disk 
in radial equilibrium. Similar enhancement is observed in the damping rate of inclinations. Numerical results are confirmed by more 
general analytical calculations that do not assume anything about origin of the radial flow in the disk. We present a simple analytical 
argument proving that the migration rate of a particle in such disk is enhanced due to the enhanced mass flux of gas colliding with 
the particle. Thus the enhancement factor does not depend on the sign of the radial gas velocity, and the migration is always directed 
inward. Within the framework of the perturbation theory we derive more general, approximate formulae for short-term variations of 
the particle semi-major axis, eccentricity and inclination in a disk out of radial equilibrium. The basic version of the formulae applies 
to the axisymmetric disk, but we present how to account for departures from axial symmetry by introducing effective components of 
the gas velocity. Comparison with numerical results proves that our formulae are correct within several percent. We have also found in 
numerical simulations that the tidal waves introduce coherence in periastron longitude and eccentricity for particles on neighbouring 
orbits. The degree of the coherence depends on the particle size and on the distance from the primary star, being most prominent for 
particles with 10 m radius. The results are important mainly in the context of planetesimal formation and, to a lesser degree, during 
the early planetesimal accretion stage. 
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1. Introduction 
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Given the number of currently known extrasolar planetary sys- 
tems (~200), it is tempting to draw some conclusions about their 
formation from statistical correlations. One of the basic correla- 
tions is the membership in a multiple stellar system. About 15% 
of the planets were found in multiple stellar systems. Assuming 
that the planets are equally common around single stars and 
in multiple systems, this is a clear discrepancy with the ob- 
served frequency of the multiplicity of solar type stars in the 
solar neighborhood, which is about 60% dDuquennov & Mayor! 
Il99ll) . A natural explanation is observational selection effects, 
because multiple systems have been usually excluded from ra- 
dial velocity searches for planets. Whether it is the only source 
of discrepancy will remain unknown until more systematic 
searches for planets in multiple syste ms are carried out. Suc h 
searches are currently underway (e.g. lMuterspaugh et al.l l2006). 
It seems reasonable, however, to expect that the planet formation 
process is influenced by companion stars. Distant companions 
(say, with periastra farther than 100 AU) do not affect planetary 
systems too much, but several extrasolar planets in binary sys- 
tems with relatively s mall separation (less tha n 100AU) have 
also been discovered (Eggenberg er et all 12004, and references 
therein). In particular, three planets have been discovered in bi- 
nary systems with separations of ~ 20 AU (HD41004A b, y 
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Cep b, G186 b). These close companions must have modified 
the structure of the protoplanetary disk (i.e. initial conditions for 
planet formation), as well as the dynamical evolution of plane- 
tary orbits. 

Although the number of extrasolar planets in multiple sys- 
tems is currently v ery small, some correlatio ns seem to be statis- 
tically significant. IZucker & Mazehl d2002l) pointed out that all 
the most massive (more than 2 Jupiter masses), short-period (pe- 
riods shorter than 40 d) planets orbit stars from binary systems. 
In consequence, for planets in binaries there is no correlation be- 
tween trms^_and_geriodaskobserved for planets around single 
stars. Eggenbe rger et al.l d2004l) showed that orbital eccentrici- 
ties tend to be very low for short-period planets. All those dif- 
ferences indicate that the companion star is affecting the planet 
formation process. 

To address planet formation in binary systems, studies of 
the orbital evolution of planetesima l s hav e been done (e.g. 
iHeppenheimeri 119781: IWhitmire et all 119981: IMarzari & Scholi 
2000; lOuintana et aT]|2002l: IThebault et a"T]|2004t IThebault et al 
2006). In the classical scenario without a gaseous disk, the 
secular perturbations from the binary companion pump up or- 
bital eccentricities of planetesimals and decelerate their accre- 
tion by reducing the gravitational focusing factor. Also, if the 
relative velocities exceed the escape velocity from their sur- 
face, collisions result in disruption rather than coagulation (e.g. 
Agnor & Asphaug 2004), s o that planetesimal accretion is in- 
hibited. Marza ri & Scholll (|2000) included the effects of uni- 
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form gas drag in an eccentric binary system and found strong 
periastron alignment of equal-size planetesimals. If the perias- 
tra are aligned, the relative velocities are kept low in spite of 
high ecce ntricities, and plane tesimal accretion is not inhibited. 
However, The bault et al.l (120061) pointed out that the alignment 
angle depends on planetesimal size, and if the size distribution 
is introduced, the relative velocities between particles of dif- 
ferent sizes are typically prohibitive for the collisional growth. 
An important result of these works is that even a small gas- 
drag force can significantly change the growth rate of planetesi- 
mals if combined with perturbations of the companion star. The 
weak point of these studies is the assumption that the gaseous 
disk is not pertu rbed, but remains stationary and axisymmetric. 
It is known (e.g.|Papaloizou & Pringlell 19771: IL"in & Papal oizoul 
1 19791: iGoldreich & T remaine 119791) . however, that the compan- 
ion induces tidal waves in the gaseous disk, which can evolve 
into strong spiral shocks. The density and velocity of the gas are 
then strongly perturbed, and the drag force acting on solid parti- 
cles is different than in the stationary, axisymmetric case. 

To date no calculations of particle motion in a disk with 
tidally induced spiral waves have been done. In a series of pa- 
pers we will explore this subject, both analytically and numeri- 
cally studying the orbital evolution and accretion of particles in 
disks perturbed by a companion star. 

In this paper we consider the case of no n- interacting particles 
orbiting in a circumprimary disk, with the perturbing companion 
star on a circular orbit. The particles range in size from 1 m to 
10 km, so our results apply to the planetesimal formation and 
early accretion phases. At first glance, the circular case may not 
seem interesting because the gravity of the companion does not 
induce secular effects on the particle orbit, like eccentricity forc- 
ing or periastra libration. However, we show that effects similar 
to those from an eccentric binary are also observed in the cir- 
cular case if the perturbations of the gaseous disk are included. 
Furthermore, we show that the orbital evolution of particles in 
such a system is significantly different than in the unperturbed, 
axisymmetric disk. The circular case is a very good starting point 
since it allows us to understand the sole effect of spiral waves in 
the gaseous disk. The eccentric binary case will be discussed in 
next paper. 

The paper is organized as follows. In Sect.[2]we present com- 
putational methods and input physics. Some important proper- 
ties of the gaseous disk are discussed in Sect. [3] In Sect. [4] we in- 
vestigate analytically and numerically the orbital evolution of a 
single particle. Section[5]is devoted to relative shapes and align- 
ment of neighbouring orbits of particles. In Sect. [6] we summa- 
rize and discuss the results. 

2. Computational method 

The problem we investigated involves the solution of both gas 
and particle equations of motion. One approach would be to 
combine hydro and N-body schemes into a single numerical 
code, however, it was enough for our purpose to perform two- 
stage simulations. We exploited the fact that in the circular bi- 
nary system the pattern of spiral waves in the gaseous disk is 
quasi-stationary in the frame co-rotating with the secondary star. 
By quasi-stationarity we mean here that the time scale of its evo- 
lution is longer than the time scale of the evolution of N-body 
particles. In the first step we obtained such a quasi-stationary 
model of the gaseous disk and then we fed it to the N-body code. 
In this way we eliminated the temporal evolution of the gas. This 
was very desirable since our goal was to investigate generic ef- 
fects of the spiral shocks on the motion of the solid bodies and 



not to perform realistic simulations of the particle growth in the 
binary system, a task we leave for future work. 

The binary system we simulated consists of primary and sec- 
ondary stars of equal mass, M p = M s = IM Q on a fixed circular 
orbit with the semi-major axis a = 23.4 AU. The implied orbital 
period is close to 80 yr. The gaseous disk and particles are or- 
biting a primary star. We chose these parameters because, apart 
from the ecce ntricity, they are c l ose to the a Centauri system 
investigated in Marzari & Scholll d2000l) . and it will enable us to 
compare their results with ours, especially in the follow-up paper 
about the eccentric binary case. Note that the self-gravity is not 
included in either fluid or particle simulations , so it is possible 
to scale the models to any size. 

2.1. Hydrodynamical simulation 

Our model of the gaseous disk was evolved using an adap- 
tive mesh refinement (AMR) code, FLASH dFrvxell et al.l 
2000 ). As a hydro solver w e employed direct Eulerian PPM 
dColella & Woodward! [T984) scheme modified to conserve an- 
gular momentum. The PPM scheme combines high-order spa- 
tial interpolation with a Riemann solver and shock-capturing 
method that results in low numerical viscosity and sharp shock 
profiles. This makes PPM particularly useful for all applications 
requiring accurate transport of momentum in a supersonic flow 
environment. 

The code solved Euler equations in 2D polar coordinates 
with the origin located at the primary star: 

— +V-(EF) = (1) 

d(YV) I r r - r r \ 

-i_^+V-(EV®V)+V/> = -Gl.\M p - + M s - + MA) (2) 

dt \ r 3 \r - r s f rj / 

where 2, P, and V = (V r , V^) denote surface density, surface 
pressure, and gas velocity at position r. The secondary star was 
located at r s such that r s = a. The equation of state was locally 
isothermal (temperature was a fixed function of distance from 
the primary star), so there was no need to solve the energy equa- 
tion, and a faster isothermal version of the Riemann solver was 
used. 

The final model of the gaseous disk was obtained as follows: 
(1) The initial disk was set up as in the case of a single star, and 
it was truncated exponentially beyond a radius slightly larger 
than the expected tidal truncation radius. (2) During the first two 
orbital periods of the binary system the grid resolution was in- 
creased up to 2048x768 in r and <p respectively (we used the 
AMR option here, but to avoid any artifacts at the edges of the 
refined blocks the whole grid was always refined). (3) The disk 
was evolved for an additional orbital period of the binary. 

The Courant number was 0.5. The grid extended radially 
from 0.4 AU to 9 AU - far enough to avoid any influence from 
boundaries on the most interesting region of the outer disk where 
the spiral waves are strongest. We took special care to minimize 
reflections at the inner and outer grid boundaries. For this pur- 
pose we tuned the standard outflow boundary conditions in order 
to reduce any discontinuities in radial direction. We also intro- 
duced a low-density hole in the first 5 radial cells at the inner 
boundary, which served as an additional buffer to damp the prop- 
agating waves. We stress that all these measures are necessary in 
order to recover the proper mass transport in the non-viscid disk. 
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2.2. Orbital integration 

Our o rbital integration code is based on Nbody4 (Aarset hTet al] 
119791) . It implements the direct summation method for self- 
gravit y, 4th-order Hermite scheme, and block time step (Makino 
1991). In this paper we consider motion of particles only un- 
der gravitational forces of the two stars and gas drag force, and 
mutual gravitational forces of planetesimals are neglected. Inter- 
particle forces and collisional accretion will be included in future 
papers. 

In the reference frame located at the primary star, the corre- 
sponding equation of motion for particle i with mass m, reads: 



The stellar configuration can simply be either a single star or 
a circular binary system. Sometimes we refer to the models 
using abbreviations presented in the Table Q] 



Table 1. Abbreviations of presented models. 





single star 


binary system 


equilibrium axisymmetric disk 


EA1 


EA2 


non-equilibrium axisymmetric disk 


NA1 


NA2 


non-axisymmetric disk 


Wl 


W2 



-G(Mp + mi)-±-GM s - 



-GM s ^+f, 



drctgji 



(3) 



where r s denotes the secondary star position. Components on the 
right hand side of Eq.[3]represent the gravity of the primary and 
secondary, the indirect term accounting for acceleration of the 
primary relative to the center of mass and gas drag force per unit 
mass, respectively. For the latter we adopt a simple formula: 



We note here that only models EA1 and W2 are physically 
consistent, while only model EA2 has been investigated by other 
authors. Thus we concentrate on differences between models 
EA2 and W2. The other models are used only as a support in 
understanding the observed effects. 



fdmgj = ~Ap\u\u, 



A = ^-C D 7Tsf, 

2m i 



3. Gaseous disk 

(4) 3.1. Disk parameters and structure 



where s, is the particle radius, p - gas density, Cd - drag coeffi- 
cient, and u is the particle velocity relative to the gas. The factor 
A is constant in our simulations. Denoting the particle's velocity 
with v we have 



u r = v r - V r , M = V - V<j. 



(5) 



We used values of Cd = 1.4 and internal density of particles 
p p = 2g/cm 3 . The Hermite scheme also requires the value of 
fdrag' which is a minor correction that was accounted for using 
numerically calculated values of p and u. 

We took special care to include gas-drag effects accurately 
yet efficiently in calculations. At the beginning of the orbital cal- 
culation, the grid data containing gas density and velocity were 
read in. During the simulation the data were rotated to match the 
current position angle of the secondary star (spiral pattern of the 
gas co-rotates with the companion), and the bi-linear interpola- 
tion was used to find the gas state at an arbitrary position of the 
particle. 

The time step was variable but limited to a maximum of 
1/(2tt • 64) yr. We tested the code with basic problems like con- 
servation of the Jacobi constant, and properly recov ered more , , 

complex results like runaw ay growth ([Kokubo & Idal [T996) and ) lr = _ - 0.05 • ( ) 

gas drag in a uniform disk dlnaba et al]|200ll) . r ' 1 AU / 



Our numerical method requires the gaseous disk to be in a state 
close to stationarity. It should also be minimally biased by nu- 
merical effects and should adequately recover all deviations from 
the Keplerian flow. To that end we used a fairly simple isother- 
mal model, which is noneth eless is close to the minimum mass 
solar nebula (lHayashilll98ll) . 
The equation of state, 



P = Sc s 2 , (6) 

is locally isothermal, with the local sound speed c s , given by the 
vertical hydrostatic equilibrium condition 

c s = -v k , (7) 
r 

where = yjGM p /r is the Keple rian v elocity and h the local 
half-thickness of the disk. Blondin(2000) has shown that, in suf- 
ficiently cold disks, the spiral waves at the outer edge of the disk 
may become unsteady. We found experimentally that for the fol- 
lowing height profile 



(8) 



2.3. The models 

Each model here is composed of the two components: gaseous 
disk and stellar system configuration. We used the following 
three gaseous disk configurations: 

- Axisymmetric disk in radial equilibrium or in a short equilib- 
rium axisymmetric disk. Characterized by a simple power- 
law density and temperature radial profiles, the radial gas 
velocity is zero. 

- Axisymmetric disk in radial non-equilibrium or in a short 
non-equilibrium axisymmetric disk. It is similar to the equi- 
librium axisymmetric disk, but the radial gas velocity is arti- 
ficially set to a non-zero value. 

- Non-axisymmetric disk, resulting from the evolution of an 
initially equilibrium axisymmetric disk in a circular binary 
system. It develops a spiral wave pattern and is naturally in 
radial non-equilibrium. 



the spiral pattern stays stable everywhere. Note that, although 
our equation of state is locally isothermal, such h r results in 
global isothermality with c s = 0.05 y/GM p /(l AU). The corre- 
sponding Mach number (vk/c s ) is 31.4 at the inner disk edge 
0=0.9 AU), which falls to 7.6 at the outer edge (r=9 AU). 

The initial profile of the surface density was given by the 
power law 



°(iiu) 



■1.5 



(9) 



which is close to the MMSN model. Since self-gravity is not 
included in hydro simulation, the normalization So can be ar- 
bitrary. To calculate the drag force (Eq. |4]i during orbital inte- 
gration, the evolved surface density £(r) was converted to the 
volume density and normalized as follows: 



p = 2 ■ 10" 



h(r)/l AU 



= 2- 10" 



'-(—] 

; UAu/ 



-1.5 



[g/cm 3 ]. (10) 
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The initial angular velocity was set to the equilibrium one for the 
given pressure gradient 



Vk + - 



1 d\nPc, 2 



2 d\nr v k 



(11) 



and the radial velocity was set to 0. 

The initial conditions described above represent our equilib- 
rium axisymmetric disk configuration. This configuration was 
evolved numerically for three binary orbital periods, which is 
enough to develop a quasi-stable spiral pattern. This evolved 
disk will be referred to as the non-axisymmetric disk. 

Here we have to point out that the 2D approximation intro- 
duces a certain inconsistency into both models of the gaseous 
disk. The problem with 2D hydrodynamical simulations is that 
the velocity field cannot be directly linked to the 3D one. The 2D 
velocity given by Eq. [TT] assures radial equilibrium for a given 
gradient of the vertically averaged pressure. Note, however, that 
it is neither the vertically averaged velocity nor velocity in the 
equatorial plane. Proper, equatorial plane velocity is given by 



1 <91n p c s 

v 'k H 

2 oln r Vk 



(12) 



where p - pc s 2 denotes pressure in the equatorial plane. Both 
formulae give different results since, in general, the gradient of p 
is different from the gradient of P. Unfortunately, the simulated 
2D velocity cannot be transformed to the equatorial plane ve- 
locity. Thus an inconsistency arises: we use 2D velocities, while 
the density is converted to 3D one (Eq.fTOli. We decided that it is 
better to use Eq.[TT]for the axisymmetric disk model although in 
principle Eq. [12] should be used: our results may not be accurate 
quantitatively, but at least we can compare both models qualita- 
tively. Furthermore we expect relative results from both models 
to be less affected than absolute ones. 

For the later considerations it is useful to define two dimen- 
sionless velocities which describe deviations from the Keplerian 
flow: 



T] = (Vk - V^)/Vk, 
K = -V r /V k . 



(13) 
(14) 



In other words, these are velocity components of a large particle 
moving on a circular, Keplerian orbit relative to the gas. 

To reveal the influence of spiral waves on the particle mo- 
tion, we have to compare results obtained in the axisymmetric 
disk with those obtained in the perturbed disk at a radius where 
angle-averaged values of p, V r , and V$ are comparable. The up- 
per panel of Fig. [TJ shows the radial profiles of tj averaged over 
the full angle for three disk models: non-axisymmetric, axisym- 
metric employing formulaQT] and axisymmetric employing for- 
mula[12] As we see, the profiles of the two first models are com- 
parable up to a distance of 3 AU from the primary star. Outside of 
this region the deviations from the Keplerian flow grow substan- 
tially. Since the spiral waves are strongest in the outer region of 
the disk, in the next section we will trace the motion of the par- 
ticle placed initially at 3 AU. The dash-dotted curve illustrates 
why the "3D velocity" is not suitable for our comparison. The 
density in the non-axisymmetric model at 3 AU has grown dur- 
ing the simulation time by roughly 15% with respect to the initial 
model (see lower panel of Fig. []]). Since the migration speed due 
to gas drag scales linearly with the density, this difference can be 
easily accounted for in comparisons with the axisymmetric disk. 




r[AU] 



Fig. 1. Top panel: radial profiles of angle-averaged r\ in non- 
axisymmetric disk (solid), axisymmetric disk with "2D veloc- 
ity" (Eq. [TT] dotted), and axisymmetric with "3D velocity" (Eq. 
[T2T dot-dashed). Bottom panel: radial profile of angle-averaged 
surface overdensity, (S - E,)/E,, in a non-axisymmetric disk. 




Fig. 2. Angular cross-sections of velocity components and den- 
sity at r-3 AU. 



3.2. Radial transport of gas 

When the spiral density waves are excited in the disk, it is no 
longer in radial equilibrium. Because the spiral pattern is rotat- 
ing slower than the local Keplerian velocity, the dissipation at 
the shock leads to a decrease in angular momentum of the or- 
biting gas and its radial mass transport. Figure [2] shows the an- 
gular cross sections of velocity components and density at 3 AU 
in the non-axisymmetric disk. Indeed, the angle-integrated mass 
flux calculated from those profiles is negative (inward). It can 
be interpreted as the result of an effective viscosity that we de- 
fine as the turbulent viscosity necessary to cause the same radial 
mass flux. In the standard a-disk theory, the kinematic viscos- 
ity v is expressed in terms of a dimensionless parameter a as 
v = 2/3ac s h, Assuming a steady accretion disk, i.e. m = -3nv"L 
independent of r, we can parametrize the mass accretion rate 
with the effective value of a e ff: 



m = 2nrV,Z — —2jra e f(C s hl.. 



(15) 



For the non-axisymmetric disk, the above equations must be 
averaged over the azimuthal angle, and finally the effective a- 
viscosity is defined as 



ffeff 



2v k V, 



eff 



3c s 



(16) 
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Fig. 3. Radial profile of the effective a-viscosity. 



where 



V, 



eff 



(17) 



It has been shown analytically dSpruitl 19871) and numerically 
dBlondinl2000l;lR6zvczka & Spruitll993b that the mass transport 
by tidal waves can be very effective. For reasonable disk parame- 
ters in close binary systems, the effective a-viscosity in the outer 
parts of the disk can easily reach 0. 1 or even more. It is hard to 
produce such high values by ordinary turbulent viscosity. 

The radial profile of a e g in our non-axisymmetric disk is 
shown in Fig. [3] Close to the inner edge of the disk, it oscil- 
lates strongly because the spiral waves are very tightly wound, 
and the radial grid resolution is insufficient for resolving them. 
However in the region close to 3 AU, the resolution is sufficient, 
and the effective a may be easily found. 



4. Evolution of orbital elements 

4.1. Non-equilibrium axisymmetric disk: analytical 
calculations 

The rate of change of orbital elements of a parti cle experiencing 
gas drag was calculated by Adac hTet al.l d!976l hereafter A76). 
Their results cannot be applied, however, applied to a disk in the 
binary system because the authors assumed that the disk is ax- 
isymmetric and stays in radial equilibrium (V r - 0). In the binary 
system the disk is neither axisymmetric nor in radial equilibrium, 
so the evolution of orbital elements may be very different. In this 
subsection, we calculate analytically the evolution of orbital el- 
ements of a particle moving through the gaseous disk, which is 
not in radial equilibrium, while still neglecting the action of the 
companion on particles. In other words, we extend the A76 work 
to the case of non-zero radial gas velocity. 

First, we consider the simplest case of the particle on a nearly 

circular, non-inclined orbit. Let u — ^juj + u^ be the value of the 

total relative velocity between particle and gas, where (u n u$) are 
radial and angular components, respectively. The particle loses 
specific angular momentum only due to angular component of 



the drag force, and for small drag (when the orbit stays nearly 
circular) its loss rate can be approximated as 



d(v^a) 1 da 



At 



2 df 



where \\ is the Keplerian velocity at radius a 
ping timescale: 



(18) 

and r is the stop- 



Apu$u 

Thus the orbit decay rate is given by the approximate formula 



(19) 



da 

— S3 -lApa — u. 
dt v k 



(20) 



There are two important factors here: and pu. The particle 
loses angular momentum when colliding with the gas at relative 
velocity u#, but the mass flux of this gas is pu, which is how the 
radial gas velocity enhances the migration rate. For large enough 
particles moving with Keplerian velocity, we can further write 



da 
~dl 



-2Apar/u. 



(21) 



In appendix A we calculate the evolution of orbital elements 
for the general case of eccentric and inclined orbits within the 
framework of perturbation theory. In the limit of a circular, non- 
inclined orbit, the general formula for orbitally averaged da/dt 
(Eq. IA.8l ) reduces exactly to Eq.l2T1 

We would like to turn the reader's attention to the depen- 
dence on the value of total relative velocity in formula ( f20b (for 
details, see Appendix). In particular, in the limit of high radial 
velocity, u r » ua, we have a very surprising relation: 



da 
~dt 



OC - M 



(22) 



The migration rate is proportional to the absolute value of gas ra- 
dial velocity; i.e., the particle migrates inward regardless of the 
direction of the radial gas flow! In order to test this result numer- 
ically we measured particle migration rates in the axisymmetric 
disk, while artificially varying the radial gas velocity. The parti- 
cle was initially on a circular orbit with components of relative 
velocity u^ and u r = V r = n ■ Us, where n was an integer num- 
ber between -5 and 5. Figure [4] shows measured and predicted 
migration rates (normalized to d(n — 1) = 1) as a function of n. 
The predicted curve is given by the formula 



d{n) 



\+n 2 



(23) 



Indeed, the migration rates measured in the simulation for n and 
—n are the same within 1%. Also the agreement with prediction 
(Eq. l23l is very good, especially for small n. For higher values 
of \n\, formula (l23l deviates slightly from results of the simula- 
tion because the orbit differs more and more from the assumed 
circular shape. In fact, changes in the semi-major axis and eccen- 
tricity are coupled so must be considered together. Even though 
the radial gas velocity does not change the angular momentum of 
the particle directly, it does change its eccentricity, which in turn 
affects the decay rate of the semimajor axis (see Eqs. 12511261 ). 

The above result concerns an idealized axisymmetric disk 
and aims to enhance basic effects of radial gas flow (regardless 
of its source), as well as to provide a proof of code validity. 
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Fig. 4. Dependence of the particle migration rate (in units of 
a(n = 1)) on the radial gas velocity. Solid line - simulation; dot- 
ted line - analytical approximation (Eq. 




axisymmetric. single 

axisymmetric, binary 

non-axisymmetric, binary 

— non-axisymmetric, binary, v R =0 




t [yr/27t] 



Fig. 5. Semimajor axis (upper panel) and eccentricity (lower 
panel) versus time for a 10 m particle in models: W2 (solid), 
EA2 (dashed), EA1 (dotted) and W2 with radial gas velocity set 
to zero (dash-dotted). 



4.2. Non-axisymmetric disk: numerical and analytical results 

In this section we present the results of the orbital calculations 
in model W2 and compare them with models EA1, EA2 and 
Wl. We also derive an analytical approximation to the perturba- 
tive formulae lA.81IA.9| for the case of non-axisymmetric disk and 
compare it with numerical results from model W2. 

As a first example, we consider particle of radius 10 m on an 
initially non-inclined, circular orbit of radius 3 AU. We checked 
empirically that in all our models the particle of this size mi- 
grates by less than 0.2 AU within the simulation time, so that the 
parameters of the gaseous disk can be regarded as roughly con- 
stant. This enables us to compare the migration speed between 
different disk models. The upper panel of Fig. [5] shows the evo- 
lution of the particle's semimajor axis in models W2, EA2, and 
EA1. As we see, the difference in migration speed of a particle 
in the non-axisymmetric disk is substantially enhanced with re- 
spect to the other two cases, roughly by a factor of three. The 
lower panel of Fig. [5] shows the corresponding evolution of ec- 



centricity. Oscillations with the synodic period of the compan- 
ion star are clearly seen. We checked that the oscillations are not 
damped by gas drag for particles larger than 10 m. The amplitude 
of oscillations depends on the semimajor axis, and the mean ec- 
centricity is roughly the same in both models with a companion 
star (only a very small decrease is observed with decreasing a). 
For this reason the eccentricity cannot be responsible for the dif- 
ference in the migration speed. We have already shown analyti- 
cally that this relatively fast migration can be induced by the ra- 
dial gas velocity. In order to test this prediction numerically, we 
performed an orbital calculation with the non-axisymmetric disk 
in which the radial gas velocity was set artificially to zero. The 
corresponding curve in Fig. [5] clearly proves that the radial gas 
flow is the main factor responsible for the accelerated particle 
migration. Another two factors are responsible for the remain- 
ing part of the difference from the axisymmetric model. First, 
in the non-axisymmetric case the effective value of if is higher 
due to weighting by the density (see discussion in the next para- 
graph). Second, due to the dissipation in spiral waves, the mean 
density in the simulation has increased by 15% with respect to 
the axisymmetric model (see Fig . |T} . 

How well do the above numerical results agree with the an- 
alytical approximations? The comparison is straightforward for 
the equilibrium axisymmetric disk. In that case formul a |A8] re- 
duces to the original formula 4.21 from A76. We have a/To = 
0.51 [2n AU/yr] in our disk model for the 10m particle on the or- 
bit with a = 3 AU. Furthermore, rj = 0.0056, and the measured 
mean eccentricity is e — 0.006. For those parameters the ana- 
lytically predicted migration velocity is within 1 % of the value 
measured in the simulation: 4.3 ■ 10~ 5 [27r AU/yr]. This proves 
the accuracy of our orbital integration. 

Application of formula lA.8l to the non-axisymmetric disk is 
not as straightforward. We have found that simply inserting an- 
gular averages of rj and k leads to quite a large discrepancy with 
the value measured in simulation. This is because the density in 
our non-axisymmetric disk is correlated with the velocity so the 
approximation IA.6I is not justified. Here we derive a new ver- 
sion of the formulae IA.8BA.9I that takes such a correlation into 
account. To that end we modify approximation |A6] by also de- 
taching p from F and averaging its product with u separately: 



(puF) = {((pu) 2 )} l/2 (F), 



(24) 



where the orbital average, () is defined by Eq. IA.5l This leads to 
the following formulae: 



t' I da 



a \dt 



r A de 
e \dt 



*eft' 



2 -2 

r + —i 

2 



1/2 



(25) 



di 
i \dt 



^ 2 I 2 , 1 -2 , 2 , 2 



1/2 



,(26) 



which are very similar to Eqs. IA.8HA.9t but variables r\ and k, 
which enter the formula for the total relative velocity (Eq. IA.4b . 
have been changed to the effective values 



2 _(pnr 

"es ~ -2 
P 



2 _(PK) 2 
*eff _ -2 ■ 
P 



(27) 



(28) 
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Fig. 6. Migration rate as a function of the particle size for mod- 
els: W2 - solid line, EA2 - dashed line, and EA1 - dotted line. 
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Fig. 7. Evolution of orbital elements for a particle on an initially 
circular, inclined orbit in models: W2 (solid line), EA2 (dotted 
line), and Wl (dashed line). 



The bar denotes here averaging over the full angle at radius r — a 
in contrast to the orbital average defined by Eq. IA.5I Also the 
characteristic time scale now depends on the averaged density: 



We note that e and i are not translated to the effective values 
because we have assumed that they remain constant during one 
orbital period. This might not be fulfilled for small particles that 
are more strongly coupled to the gas. Fortunately, the effect be- 
comes noticeable only at the lower limit of the size range for 
which the drag law we used is applicable. In analogy to Eq.l28l 
we can define the effective a-viscosity, and its relation to K e ff is 

*eir = \h 2 r A^J- (30) 

Now we are in a position to test this result with numerical 
simulations. From the model W2 we measured mean e = 0.006, 
77 = 0.0056, and effective K e ff = 0.027, rj e ff - 0.01. For these pa- 
rameters, the predicted migration velocity is 1 .9 TO -4 [27r AU/yr] 
in comparison to 1.8 • 10 -4 [27r AU/yr] measured from the simu- 
lation. Taking the number of approximations that had to be done 
in the analytical derivation into account, this agreement is indeed 
excellent. 

In Fig. |6]we display the migration rates measured in simula- 
tions as a function of the particle size. As we see, the presence 
of the companion star does not change the migration rate in ax- 
isymmetric disks, because the excited eccentricity oscillations 
are small in comparison to the 77 parameter that plays the major 
role. In the non-axisymmetric disk parameter /e e jf clearly dom- 
inates other factors (e, i, 77^), and thus for all considered sizes 
migration speed is enhanced by roughly the same factor of 3. 
Of course this factor will change as K e g is changing. The value 
of /c e ff depends on the position in the disk, disk properties, and 
binary separation. In particular, for a given disk model, it de- 
creases with decreasing r, because the waves are weaker in the 
inner disk. For a discussion of the dependence of a e ff on the disk 
density profile a nd temperature (for isothermal models), please 
refer to lBlondinld2000l) . 

To test the formula for inclination damping (Eq. |26| i, we 
placed a particle on an initially circular orbit at 3 AU with in- 
clination 0.01. The evolution of its orbital elements in models 



Wl, W2, and EA2 is shown in Fig. [7] Clearly, the inclination 
decreases faster in the non-axisymmetric disk. The initial differ- 
ence in damping rate, measured over interval 5t = 5Qyr/2n, is 
about three times higher than in model EA2. In order to make 
a comparison with the analytical approximation (which does not 
account for the presence of a secondary), we ran one more model 
with the non-axisymmetric disk but with the gravity of the sec- 
ondary star switched off (Wl; see Fig. [5). The difference be- 
tween models Wl and W2 is very small, which is expected since 
the companion star influences the inclination rather weakly. We 
measured the initial slope in model Wl to be 2.2 • 10~ 5 [2n/yr], 
while the approximate Eq. [26] gives 2.5 ■ 10~ 5 [2n/yr], which 
makes a 12% difference. Althrough worse than for the semi- 
major axis, this accuracy is still very good for an approximate 
formula. 

The eccentricity of the particle in a binary system is set by 
the balance between forcing from the secondary and drag damp- 
ing. As it was already shown in Figure [5] the eccentricity of 
a particle on an initially circular orbit is immediately forced 
to oscillate with a roughly constant amplitude (0.01 at 3 AU). 
Comparison of the analytical formula for damping of e (Eq. |26| | 
with simulation results in model W2 would require setting up the 
particle with initial eccentricity substantially higher than 0.01. 
Since it is rather unlikely for a small particle to have such eccen- 
tricities, and in addition analytical approximation may not work 
well in this regime, we will not present it. We note only that the 
comparison of eccentricity between models Wl and W2 in Fig. 
Qreveals that the spiral waves are only responsible for the lower 
limit of eccentricity in model W2, while the remaining contribu- 
tion comes from the perturbation by the secondary. 

5. Coherence of periastra and eccentricities 

Relative shape and alignment of neighbouring orbits is a very 
important factor because it controls the relative velocities of the 
particles and thus their growth rate. Even if the companion star 
excites high eccentricities, it does not necessarily m ean that the 
relative velocities are high.|M arzan & Scholll fcOOO) have shown 
that secular perturbations from an eccentric companion, com- 
bined with the gas drag, lead to the strong alignment of periastra 
between particles of the same size on neighbouring orbits. Since 
periastra are also coupled to eccentricities, the relative veloci- 
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Fig. 8. Periastron longitude a> in model W2. Each row shows Fig 9 Eccen tricity as a function of semimajor axis. Each row 
temporal evolution for a particle of a given size. The time is dis- shows the temporal evolution for a particle of a given size, 
played in years. 



ties are low. This effect of "orbital phasing" in an eccentric bi- 
nary is prominent even for bodies of 100 km in diameter, which 
are commonly regarded as decoupled from the gas. It should be 
stressed that the collisio n velocities are low only between parti- 
cles of the same size. [Thebault et al. (2006) has shown that for 
distribution of sizes, even small misalignment in the periastra of 
particles of different sizes results in relative velocities that are 
high enough to prohibit collisional growth (note however that 
they used an equilibrium axisymmetric-disk approximation). 

In this section we investigate the effects of the orbital coher- 
ence, but in the circular binary system and for smaller sizes of 
particles. The secular effects of eccentricity forcing and periastra 
alignment are not present in the circular system. However, the in- 
clusion of a non-axisymmetric gaseous wave pattern provides an 
additional factor that localy perturbs orbit of the particle. Since 
the wave pattern co-rotates with the companion star, we may 
expect similar effects to the secular gravitational perturbations. 
To check what really happens in such systems, we performed 
4 runs with non-interacting particles of sizes 1 m, 10 m, 100 m, 
and 1 km. Each run was initiated with 30000 particles on cir- 
cular, non-inclined orbits distributed randomly between 0.8 AU 
and6AU. 

Figure [8] shows the longitudes of particle periastra with re- 
spect to the longitude of the companion, u>, as a function of par- 
ticle semimajor axis. Each row shows the evolutionary sequence 
for the population of particles of different size (indicated on the 
vertical axis legend). The time evolution is shown to prove that 
the observed configurations are not transient but converge to a 
certain, stationary pattern on the a> - a plane. 

A quick look at the plot reveals that the effect of periastra 
alignment is present for all considered sizes of particles. The 
degree of alignment depends, however, on the distance from the 
central star and on the particle size. A closer exploration of the 
plot allows two different regimes to be distinguished depending 
on the particle size: 

1 . Particles with radii smaller than a few meters feel strong gas 
drag, and their periastra are correlated with the spiral pattern 
in the gaseous disk. It can be observed on the plot for 1 m 
particles in the region of stable orbits (below 3 AU). 



2. Particles with radii larger than a few meters have orbits 
aligned, but the alignment is not correlated with the waves 
in the gaseous disk. In the outer parts of the disk, around 
3 AU, the periastra are in opposition to the companion star 
with a substantial scatter in longitude (larger for larger par- 
ticles). When the drag force is stronger (in the inner disk or 
for smaller particles), the alignment is more pronounced. It 
is particularly strong for 10 m particles. 

We want to stress here that the alignment of periastra does 
not imply any spatial alignment of particles. In fact we did not 
observe a spatial correlation of particles with the spiral waves 
at the particle-size range considered here. We only checked that 
such a correlation becomes weakly visible for 10 cm particles. 

Figure [9] illustrates the dependence of eccentricity on semi- 
major axis in the same manner as for periastra longitudes. The 
coherence in eccen tricity is much weaker in comparison to the 
eccentric model of Mar zari & Scholll d2000j) . In the inner disk, 
it is simply the effect of strong damping by the gas drag. The 
presence of spiral waves manifests mainly for 1 m and 10 m size 
particles in the form of pulsations in a. Actually, two pulsation 
patters can be noticed for 1 m particles corresponding to two spi- 
ral arms. We expect that the relative phase between those two 
patterns depends on the winding angle of the spirals. Thus the 
coherence in eccentricity will vary depending on the disk tem- 
perature. 

6. Discussion and summary 

We have investigated orbital evolution of particles moving in a 
circumprimary gaseous disk in a circular binary system. This is 
the first investigation of the orbital motions that includes tidally 
excited density waves in the gas-drag calculation. 

We have demonstrated numerically that the radial flow of 
gas in the disk (inward or outward) effectively increases the drag 
force due to an enhanced mass flux colliding with the solid parti- 
cles. We derived approximate analytical formulae for the change 
rate of orbital elements (a, e, i) in the gaseous disk, which is 
not in radial equilibrium. The formulae do not assume anything 
about the source of the radial gas flow. In general there are four 
constituents of the relative velocity between particle and gas. 
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The first two, eccentricity and inclination, describe the devia- 
tion of the particle from the Keplerian gas flow. The remaining 
two represent the radial and angular deviation of the gas from the 
Keplerian flow. To account for the non-axisymmetric features in 
the disk, one has to introduce effective components of the gas ve- 
locity. The effective components can be substantially larger than 
simple angular averages, meaning that the evolution of orbital el- 
ements can be faster in the non-axisymmetric disk. Exactly this 
situation happens in the investigated disk perturbed by the com- 
panion star: the tidally induced spiral waves which propagate 
radially in the form of shocks increase the effective components 
of the gas velocity. In particular, in the outer disk where the per- 
turbation is strongest, the effective radial component dominates 
other components of the relative velocity. We have found nu- 
merically that the particles of sizes from 1 m to 10 km migrate 
there around three times faster than in an axisymmetric disk in 
radial equilibrium. This result agrees with the analytical predic- 
tion very well, within a few percent. 

The effect of enhanced drag force due to radial gas flow 
has some significant consequences in the context of planetesi- 
mal formation in the binary system. The faster particle migra- 
tion raises the old question of whether there is enough time to 
form planetesimals before all smaller particles have fallen onto 
the star. We have to note, however, that the enhancement in the 
drag force gradually vanishes inward in the disk. Depending on 
the disk temperature and density radial profile, this differential 
migration may eventually result in the accumulation of particles 
at a certain radius in the inner disk, although we estimate that 
this is not the case for realistic disk models. There are other ef- 
fects that may support faster planetesimal formation. If only the 
larger bodies can form quickly, the size-independent enhance- 
ment of the migration speed will result in a higher flux of smaller 
particles, which will feed the larger body. Furthermore, an accel- 
erated damping of the particle's inclinations leads to an increase 
in their number density in the mid-plane of the disk and more 
frequent collisions. Both effects are important for the planetesi- 
mal formation and early stages of their growth. 

The frequency of collisions between particles is only one of 
the two factors controlling their growth rate. The other one is the 
relative velocity of colliding bodies that determines whether the 
bodies stick together or shatter into smaller pieces. The relative 
velocity almost exclusively depends on the eccentricity and pe- 
riastra orientation of the particles, provided that the inclinations 
are already damped. The contribution from decay of the semi- 
major axis is negligible (in our system it can play a role only 
for 1 m particles at 3 AU). We find that the spiral waves induce a 
certain coherency in both periastra longitudes and eccentricities. 
Only for the smallest particles of 1 m size does the coherence in 
periastra longitudes come from a simple correlation of the or- 
bit orientation with the spiral wave pattern in the gas. For larger 
bodies, there is no such correlation, and the degree of coher- 
ence decreases with increasing size of the body. It seems that the 
relative velocities will be affected by coherence only for bodies 
smaller than 10 m. In order to determine how much it does help 
in planetesimal formation, collisional simulations with realistic 
size distribution will have to be carried. This will be the subject 
of the next paper. 

In this paper we focused on the circular binary system be- 
cause it allows the sole effects of the spiral waves to be studied 
without time dependent effects that are present in the eccentric 
binary. On the other hand, the radial gas flow in the eccentric 
system is much stronger and we expect it to influence particle 
motion to a much higher degree than in the circular case. The 



eccentric binary case will be investigated in one of the subse- 
quent papers. 
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Appendix A: Appendix: Mean variations of orbital 
elements 

In this section we extend calculations of A76 for the case of a 
disk in radial non-equilibrium. All other assumptions made in 
A76 are preserved; in particular, the disk is axisymmetric. Parts 
of the formulae written in bold font are additions with respect to 
their formulae. 

Introducing non-dimensional radial velocity k = — V r /vk, 
Eqs. (4.7) in A76 take the following form: 



da 
dt 

de 
dt 

di 
dt 



-Apu (l + 2e cos if/ + e 2 - (1 + e cos if/) 3/2 h cos i + 

1 - e 1 1 

+ *■(!- e 2 ) 1/2 e sin if/] , (A.l) 



2 cos ip + e + e cos 2 ip 
-Apu < 2 cos ip + 2e — — — h cos i+ 



(1 + e cos if/) 



1/2 



+ 41 - e 2 ) 1 ' 2 sin if/] , 
h 



= -Apu 



(1 + e cost/0 



1/2 



cos (ip + to) sin i, 



(A.l) 
(A.3) 



where if/ and u denote true anomaly and argument of periastron 
respectively, and h » 1 — 77. The total relative velocity u reads: 

^000 00 

u — Vk (a){(l — — cos if/)e + cos (ip + a>)i + r\ + 77e cos t/>} + 



l-e : 



{**(! - e 2 ) - 2/ce(l - e 2 ) 1 ' 2 sin if/] . (A.4) 



Assuming that the drag force is weak and the orbital elements 
are constant during one orbital period, their orbitally averaged 
rate of variation is expressed as 



dQ 
dt 



1 

2n j l} 



2k 



dQ (l-e 2 ) 3/2 
dt (1 + ecosi^) 1 



dip, 



(A.5) 



where Q e {a, e, 2'}. In order to simplify further calculations we 
apply the approximation (Eq. 4. 19 in A76): 



(uF) = {u)(F) = [(u 2 )) XI2 (F) 



(A.6) 



where F denotes factors other than u in Eqs. ( IA. lt -( IA.3l >. It is 
justified as long as u and F are independent and the constant 
part of u 2 is greater than the amplitude of the oscillatory part. 
Providing that variables e, i, 77, k are much smaller than unity and 
preserving only leading terms for each of them, we obtain from 
Eq. (E3) 



(u 2 ) = v k 2 (a) 



1 



-i 2 + i 1 2 + l ? 



(A.l) 



In practice, term k 2 in front of e 2 is negligible and the remain- 
ing formula is intuitive: in addition to the angular component of 
velocity, 77, its radial part, k, enters. 
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Fortunately, the orbital average of F is not changed with 
respect to A76, since all new terms in Eqs. ( [A.ll i- dATJt vanish 
when averaged. Thus the only changes come from < u 2 > term 
and finally we obtain: 



a \dt 



9 1 .9 1 ? 

- -iC\e + -i +rj + k 



1/2 



(A.8) 



to /de\ _ 2 l2.l—\ - 
e \ df / i \ dt 



— — K\e +7 ? + K 



1/2 



(A.9) 
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